Datos observacionales y estrategias de identificación
Ayudantía 2- Curso métodos en CCS
I. Bienvenida al tutorial
Hola, este tutorial está escrito como la actividad de aplicación de la ayudantía 2 del curso
La idea es ir desarrollando el tutorial en paralelo a la sesión y que puedan tomar apuntes adicionales directamente en el notebook. Está escrito como .qmd, un notebook Quarto y puede ser editado en rstudio, vscode o cualquier otro IDE o editor de texto. El lenguaje utilizado para el procesamiento y análisis es R.
Se recomienda crear un fork al repositorio y trabajar en su propia rama.
Como funciona
El tutorial se estructura como un proyecto de investigación ficticio completo, que empezaremos cargando datos brutos hasta obtener e interpretar coficientes de regresión.
Lo métodos a cubir son:
1. Regresión lineal y método de diferencias (en datos experimentales)
[3. Variables Instrumentales]
[4. Regresión Discontinua]
[5. Diferencias en Diferencias]
[6. Efectos fijos, aleatorios y modelos jerárquicos]
[7. Matching]
Los datos y notebook son descargables, así que se recomienda correr todo en tu propio computador o en una instancia en la nube.
Este tutorial es una adaptación del trabajo de Hans H. Sievertsen, con la adición de métodos de matching y con cambios a la situación a analizar.
1 Pregunta de investigación y datos
Vamos a usar una situación ficticia, con datos simulados para poder aplicar de manera consisa los diferentes tipos de estrategias de identificación revisadas en clase. ALgunas (que se supone ya manejan) las revisaremos muy rápidamente y en otras, vamos a tener mayor énfasis.
1.1 Pregunta de investigación
Nuestro objetivo es responder la siguiente pregunta ficticia de investigación:
Asistir a preuniversitario mejora los resultados académicos?
Para responder esta pregunta, usaremos unos datos ficticios y simulados
1.2 Contexto ficticio
La pregunta de investigación se inspira en trabajos como el de Matsudaira (2007) e intervenciones en estudiantes de bajo nivel socioeconómico por Dietrichson et al ( 2017).
El escenario ficticio es el siguiente:
- Para un conjunto de colegios en una comuna, existe la opción de asistir a pre-unirversitario intensivo los fines de semana durante el segundo semestre de 4to medio.
- El preuniversitario se enfoca en preparar la prueba de admisión a la universidad vigente (PSU en ese momento)
- El preuniversitario es gratuito, pero para ser matriculados requiere que los padres se involucren en un proceso.
- Estamos interesados en testear el impacto de la participación en el preuniversitario en los resultados académicos de los estudiantes.
1.3 Datos ficticios dispobibles
- school_data_1.csv
- Usamos esta data para ejemplificar como cargar data guardada en formato csv.
- Este dataset tiene información sobre cada individuo (con identificador id), la escuela a la que asiste, un indicador si participó en el preuniversitario, sexo, ingreso del hogar, educación de los padres, resultados en diagnostico PSU a incios del año escolar y el puntaje efectivamente obtenido en la prueba oficial.
- school_data_2.dta
- Usamos esta data para ejemplificar como cargar data guardada en formato STATA.
- Este dataset tiene información de cada individuo (con identificador id).
- Este dataset tiene la información si el individuo recibió la carta de invitación para participar del preuniversitario.
- school_data_3.xlsx
- Usamos este dataset para practicar como cargar datos guardados en formato Microsoft Excel.
- Este dataset incluye datos sobre cada individuo (con identificador id)
- Este dataset tiene información de rendimiento académico antes y después del preuniversitario.
II. Preparación de los datos
En esta sección vamos a preparar los datos para el análisis. Este proceso generalmente incluye cargarlos, inspeccionarlos, limpiar y dar la estructura deseada. También revisaremos estadísticas descriptivas que nos den una idea antes de estimar cualquier modelo.
En la vida real, este proceso suele ser bastante largo, laborioso e implica volver sobre lso pasos anteriores múltiples veces. También invlocura tomar desiciones por parte de los investigadores, por lo cual la documentación de esta fase es especialmente importante.
En nuestro caso, será bastante lineal y directo, ya que son datos ficticios simulados. Pero en la realidad, no suele ser así.
1 Cargar los datos
1.1 Intalar y cargar paquetes
Para poder cargar los datos, necesiamos los paquetes adecuados. Siempre hay multiples formas de hacer las cosas en cualquier lenguaje o programa. En este caso, usaremos la función read_csv() del paquete readr. Para poder usarlo, debemos estar seguros de que está instalado.
Si no está instalado, podemos hacerlo con la función install.packages("[nombre paquete a instalar]")
install.packages("readr")Si el paquete ya está instalado, para poder usarlo necesitamos tenerlo cargado en nuestra librería. Para esto usamos la función library("[nombre paquete a cargar]").
library(readr)Cada paquete lo debemos instalar solo una vez por computador, pero debemos cargarlo en cada sesión para poder utilizarlo.
1.2 Cargar datos csv
Con readr estamos en condiciones de usar la función read_csv() para cargar la primera base de datos.
Vamos a cargar school_data_1.csv agregando el path a los datos en paréntesis. Puede reemplazar por el path correspondiente o usar el paquete auxiliar here
PS. notemos que tambien estamos incluyendo comentarios en los bloques de código. Las líneas que empiezan con el símbolo # son ignoradas por R.
# load readr package
library("readr", "here")
#prueba el directorio que te dice con here
here::here()[1] "/Users/melita/Documents/GitHub/ayudantiasccs"Entonces, dentro de here("[Escribes el directorio relativo]") actuará como el directorio relativo sin errores. Entonces, cargamos los datos y los asignamos a un data frame.
school_data_1 <- read.csv(here::here("data_raw/school_data_1.csv"))Usualmente, después de cargar un dataset es útil visualizarlo. Empleamos la función head() para ver sus primeras 6 observaciones.
head(school_data_1) person_id school_id summercamp female parental_schooling parental_lincome
1 1 5 0 1 10 12.93921
2 2 14 1 0 11 14.73739
3 3 7 1 0 14 16.08326
4 4 8 0 0 12 14.60808
5 5 9 1 0 11 13.75679
6 6 26 1 1 11 14.70552
test_year_5 test_year_6
1 NA 1.636284
2 1.297036 2.561449
3 2.761908 3.526002
4 1.772913 2.381062
5 1.180356 1.959818
6 2.439802 2.810571Una segunda alternativa es descargar y cargar los datos directamente en R desde internet. Puede ser cualquier link directo o, si está alojado en github, tienes que asegurarte de que sea un repositorio público.
school_data_1 <- read.csv("https://github.com/ClasesMOW/ayudantiasccs/blob/main/data_raw/school_data_1.csv")1.3 Cargar datos STATA
En Ciencias Sociales y Economía es muy comun contar con datos para ser utilizados en el programa STATA, estos son archivos que terminan en .dta. Para cargarlos, vamos a usar HAVEN del paquete tidyverse. También usaremos muchas otras funciones de ese paquete, asi que es buen momento para cargalo.
library(tidyverse)
school_data_2<- haven::read_dta(here::here("data_raw/school_data_2.dta"))Usemos el comando tail() para ver las últimas 10 entradas.
# print the 8 rows in the tail
tail(school_data_2,n=10)# A tibble: 10 × 2
person_id letter
<dbl> <dbl>
1 3482 1
2 3483 0
3 3484 0
4 3485 0
5 3486 1
6 3487 0
7 3488 1
8 3489 0
9 3490 0
10 3491 01.4 Cargar datos Microsoft Excel
Finalmente, cargaremos el tercer dataset guardado como hoja de cálculo de Excel .xlsx. Para esto usaremos el paquete readxl que viene incluido en tidyverse. Luego de asignarlo a un dataframe, démosle una mirada con glimpse() (también del tidyverse).
school_data_3 <- readxl::read_xlsx(here::here("data_raw/school_data_3.xlsx"))
glimpse(school_data_3)Rows: 3,491
Columns: 10
$ person_id <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17…
$ test_year_2 <dbl> 1.090117, 1.246309, 2.726472, 2.693032, 1.660545, 2.22377…
$ test_year_3 <dbl> 1.914594, 1.154470, 2.269011, 2.413203, 1.828067, 2.27566…
$ test_year_4 <dbl> 2.065805, 1.582455, 3.247252, 1.479452, 1.361972, 2.38510…
$ test_year_7 <dbl> 2.377697, 1.747376, 3.017764, 2.637954, 1.904636, 3.37613…
$ test_year_8 <dbl> 2.032904, 2.444041, 3.361646, 3.021940, 2.109774, 3.24542…
$ test_year_9 <dbl> 1.493803, 1.663050, 3.387020, 2.761513, 2.285818, 2.96503…
$ test_year_10 <dbl> 1.880512, 1.833769, 2.968617, 2.088086, 1.845694, 3.30819…
$ learnings <dbl> 10.236394, 8.278911, 8.966529, 8.876466, 8.770518, 10.484…
$ school_id <dbl> 5, 14, 7, 8, 9, 26, 13, 11, 23, 9, 25, 15, 3, 4, 17, 7, 1…2. Unir los datasets
Tenemos 3 bases de datos con información diferente de los mismos individuos. Generalmente es buena idea tener una sola gran tabla con toda esta información, especialmente si estimaremos modelos en base a ésta.
La base de datos 1 y 2 tienen una forma similar: los individuos son filas y las variables columnas y hay una sola fila para cada individuo.
Para hacerlo, podriamos usar varias alternativas.
Una alternativa es usar la función nativa merge( ). En esta función, primero mencionamos los datasets a unir, luego informamos cual es la(s) columnas(s) que debe usar para unir ambos datasets con by=..... Por defecto, R incluye todas las filas que están en ambos datasets (basados en la variable by), pero podemos fijar all=TRUE para mantener todas las filas que están en ambos datasets o all.x=TRUE para mantener todas las filas coincidentes y las del primer dataset or all.y=TRUE para guardar todas las filas del segundo dataset.
Veamos un ejemplo con los dos primeros datasets. Luego usemos dim() para conocer las dimensiones del nuevo dataset unido en términos de filas y columnas.
# merge school_data_1 with school_data_2
school_data_merged <-merge(school_data_1,school_data_2, by="person_id")
# dimensions
dim(school_data_merged)[1] 3491 9Notemos que el dataset unido tiene 3491 filas y 9 columnas. Unimos todas las filas y agregamos al dataset de información del estudiante si recibió o no la carta (school_data_2)
¿Qué ocurre si las columnas tienen igual nombre? R va a renombrarlas automáticamente agregando un sufijo .x (a la columna del primer dataset) e .y ( a la columna del segundo dataset).
Entonces, en el siguiente bloque:
- Unimos school_data_1 y school_data_2 usando como variable de unión person_id y guardamos el dataset unido como school_data.
- Unimos school_data_3 con school_data y sobre-escribimos school_data.
Notar que acá unimos por las columnas person_id y school_id. Esto no es realmente necesario porque cada estudiante con id única tiene un solo colegio, pero sirve de ejemplo en como usar más de una columna mediante c(). 3. Usamos la función summary() para obtener una estadística descriptiva de las variables en el dataset unido.
# Merge school_data_1 and school_data_2 and save as school_data_merged
school_data_merged<-merge(school_data_1,school_data_2,by="person_id")
# Merge school_data_3 with school_data_merged
school_data_merged<-merge(school_data_merged,school_data_3,by=c("person_id","school_id"))
# summary statistics
summary(school_data_merged) person_id school_id summercamp female
Min. : 1.0 Min. : 1.00 Min. :0.0000 Min. :0.0000
1st Qu.: 873.5 1st Qu.: 8.00 1st Qu.:0.0000 1st Qu.:0.0000
Median :1746.0 Median :15.00 Median :0.0000 Median :1.0000
Mean :1746.0 Mean :15.66 Mean :0.4589 Mean :0.5153
3rd Qu.:2618.5 3rd Qu.:23.00 3rd Qu.:1.0000 3rd Qu.:1.0000
Max. :3491.0 Max. :30.00 Max. :1.0000 Max. :1.0000
parental_schooling parental_lincome test_year_5 test_year_6
Min. :10.00 Min. :12.67 Min. :-0.2696 Min. :0.2824
1st Qu.:11.00 1st Qu.:14.12 1st Qu.: 1.7997 1st Qu.:2.0071
Median :11.00 Median :14.52 Median : 2.2583 Median :2.4492
Mean :11.32 Mean :14.56 Mean : 2.2589 Mean :2.5174
3rd Qu.:12.00 3rd Qu.:14.95 3rd Qu.: 2.7154 3rd Qu.:3.0185
Max. :23.00 Max. :19.45 Max. : 4.6897 Max. :4.9831
NA's :5 NA's :6 NA's :5
letter test_year_2 test_year_3 test_year_4
Min. :0.0000 Min. :-0.5701 Min. :-0.464 Min. :0.01425
1st Qu.:0.0000 1st Qu.: 1.7712 1st Qu.: 1.818 1st Qu.:1.78282
Median :0.0000 Median : 2.2377 Median : 2.249 Median :2.24947
Mean :0.2469 Mean : 2.2246 Mean : 2.254 Mean :2.24952
3rd Qu.:0.0000 3rd Qu.: 2.6644 3rd Qu.: 2.715 3rd Qu.:2.70903
Max. :1.0000 Max. : 4.7373 Max. : 4.389 Max. :4.34995
test_year_7 test_year_8 test_year_9 test_year_10
Min. :0.03139 Min. :-0.1856 Min. :-0.2456 Min. :-0.3043
1st Qu.:1.93768 1st Qu.: 1.9277 1st Qu.: 1.9371 1st Qu.: 1.9415
Median :2.36855 Median : 2.3889 Median : 2.3868 Median : 2.3776
Mean :2.44445 Mean : 2.4506 Mean : 2.4531 Mean : 2.4541
3rd Qu.:2.97207 3rd Qu.: 2.9596 3rd Qu.: 2.9861 3rd Qu.: 2.9693
Max. :4.77111 Max. : 5.0451 Max. : 4.6565 Max. : 4.9295
learnings
Min. : 5.386
1st Qu.: 9.181
Median :10.041
Mean :10.064
3rd Qu.:10.959
Max. :14.621
Otra opción de hacer lo mismo es con join del paquete dplyr del tidyverse.
(foto sacada de link, puedes ver ejemplos detallados de los diferentes tipos de join)
En este caso, podríamos hacer en un pipe %>% dos left joins seguidos.
school_data_merged <- school_data_1 %>%
left_join(school_data_2, by="person_id") %>%
left_join(school_data_3, by=c("person_id", "school_id"))3. Lipiar los datos
3.1 Tidyng los datos
Ahora que hemos unido las bases de datos, trataremos de que satisfazgan los principios de Tidy Data.
A data frame is said to be tidy (Hadley’s paper) if the following conditions are met:
- A single column (variable) contains all measures of the same data feature of the observations.
- A row contains measures of the same observation.
(puedes profundizar y ver más ejemplos aplicados en https://sscc.wisc.edu/sscc/pubs/DWE/book/6-2-tidy-data.html )
Uno de estos, es que cada columna debe ser una variable y cafa fila una unidad de observación.
Si inspeccionamos el número de columnas:
# Alternatively we can use the nrow() functions to get the number of columns
ncol(school_data_merged)[1] 17Son 17, pero tenemos 9 variables.
head(school_data_merged) person_id school_id summercamp female parental_schooling parental_lincome
1 1 5 0 1 10 12.93921
2 2 14 1 0 11 14.73739
3 3 7 1 0 14 16.08326
4 4 8 0 0 12 14.60808
5 5 9 1 0 11 13.75679
6 6 26 1 1 11 14.70552
test_year_5 test_year_6 letter test_year_2 test_year_3 test_year_4
1 NA 1.636284 0 1.090117 1.914594 2.065805
2 1.297036 2.561449 0 1.246309 1.154470 1.582455
3 2.761908 3.526002 0 2.726472 2.269011 3.247252
4 1.772913 2.381062 0 2.693032 2.413203 1.479452
5 1.180356 1.959818 0 1.660545 1.828067 1.361972
6 2.439802 2.810571 0 2.223775 2.275665 2.385106
test_year_7 test_year_8 test_year_9 test_year_10 learnings
1 2.377697 2.032904 1.493803 1.880512 10.236394
2 1.747376 2.444041 1.663050 1.833769 8.278911
3 3.017764 3.361646 3.387020 2.968617 8.966529
4 2.637954 3.021940 2.761513 2.088086 8.876466
5 1.904636 2.109774 2.285818 1.845694 8.770518
6 3.376130 3.245426 2.965039 3.308194 10.484683Es porque tenemos puntajes de las pruebas del año 2 al 10. Este tipo de datos son de panel
Generalmente, que hagamos con este tipo de datos depende del tipo de modelos que queramos usar. Si bien el formato wide es facil de entender, generlamente para modelos y análisi preferimos que esté en formato long. Especialmente cuando modelamos incluyendo efectos fijos También es este el que adhiere a los principios tidy de mejor manera.
Para cambiar a long, usamos pivot_longer()
# load tidyr package. Ya etsá incluido en tidyverse, pero también se puede llamar por si solo.
library("tidyr")
# make data tidy (make long)
school_data_tidy<-school_data_merged%>%
tidyr::pivot_longer(
cols = starts_with("test_year"),
names_to = "year",
names_prefix = "test_year_",
names_transform = list(year = as.integer),
values_to = "test_score",
)
# ncol to get the number of columns of the new dataset
ncol(school_data_tidy)[1] 10Ahora tenemos nuestros datos listos para que los inspeccionemos.
3.2 Selección de muestra
Ya que contamos con datos que siguen los principios de tidy data, lo siguiente es seleccionar la muestra apropiada. En este trabajo, los unicos problemas que podríamos enfrentar son relacionados con valores faltantes o missing. Para inspeccionarlos vamos a usar la función skim() del paquete skimr, esta función nos muestra los vaores faltantes en nuestro dataset de una manera global.
# Cargar skimr
library("skimr")
# Usamos skim() para inspeccionar los datos
skim(school_data_tidy)| Name | school_data_tidy |
| Number of rows | 31419 |
| Number of columns | 10 |
| _______________________ | |
| Column type frequency: | |
| numeric | 10 |
| ________________________ | |
| Group variables | None |
Variable type: numeric
| skim_variable | n_missing | complete_rate | mean | sd | p0 | p25 | p50 | p75 | p100 | hist |
|---|---|---|---|---|---|---|---|---|---|---|
| person_id | 0 | 1 | 1746.00 | 1007.78 | 1.00 | 873.00 | 1746.00 | 2619.00 | 3491.00 | ▇▇▇▇▇ |
| school_id | 0 | 1 | 15.66 | 8.67 | 1.00 | 8.00 | 15.00 | 23.00 | 30.00 | ▇▇▇▇▇ |
| summercamp | 0 | 1 | 0.46 | 0.50 | 0.00 | 0.00 | 0.00 | 1.00 | 1.00 | ▇▁▁▁▇ |
| female | 0 | 1 | 0.52 | 0.50 | 0.00 | 0.00 | 1.00 | 1.00 | 1.00 | ▇▁▁▁▇ |
| parental_schooling | 45 | 1 | 11.32 | 1.10 | 10.00 | 11.00 | 11.00 | 12.00 | 23.00 | ▇▁▁▁▁ |
| parental_lincome | 0 | 1 | 14.56 | 0.69 | 12.67 | 14.11 | 14.52 | 14.95 | 19.45 | ▂▇▁▁▁ |
| letter | 0 | 1 | 0.25 | 0.43 | 0.00 | 0.00 | 0.00 | 0.00 | 1.00 | ▇▁▁▁▂ |
| learnings | 0 | 1 | 10.06 | 1.28 | 5.39 | 9.18 | 10.04 | 10.96 | 14.62 | ▁▃▇▃▁ |
| year | 0 | 1 | 6.00 | 2.58 | 2.00 | 4.00 | 6.00 | 8.00 | 10.00 | ▇▇▃▇▇ |
| test_score | 11 | 1 | 2.37 | 0.72 | -0.57 | 1.88 | 2.33 | 2.84 | 5.05 | ▁▂▇▃▁ |
Con esta función podemos ver facilmente cuantas filas y comunas son, los tipos de varables y número de missing values. Además la media, desviación estándar, percentiles e incluso un histograma para cada variable.
En estos casos, para parental_schooling tenemos 45 missing y para test_score 11. Asumamos que estos valores missing son random y deseamos remover estas filas. Para esto usamos filter(). Esta funcion toma dos argumentos, el dataset a filtrar y la condición para que se mantenga en el dataset, en este caso que no sea na o !is.na(partental_schoolin). La función is.na() es verdad cuando el elemento en () es missing y usamos !para mostrar que queremos lo contrario a esto condición. En otras palabras queremos que la educaión parental no esté missing.
# Seleccionamos las columnas sin missing values
school_data_selected<-dplyr::filter(school_data_tidy,!is.na(parental_schooling),!is.na(test_score))
# Use skim() to skim the data again
skim(school_data_selected)| Name | school_data_selected |
| Number of rows | 31363 |
| Number of columns | 10 |
| _______________________ | |
| Column type frequency: | |
| numeric | 10 |
| ________________________ | |
| Group variables | None |
Variable type: numeric
| skim_variable | n_missing | complete_rate | mean | sd | p0 | p25 | p50 | p75 | p100 | hist |
|---|---|---|---|---|---|---|---|---|---|---|
| person_id | 0 | 1 | 1746.67 | 1007.89 | 1.00 | 873.00 | 1746.00 | 2620.00 | 3491.00 | ▇▇▇▇▇ |
| school_id | 0 | 1 | 15.66 | 8.68 | 1.00 | 8.00 | 15.00 | 23.00 | 30.00 | ▇▇▇▇▇ |
| summercamp | 0 | 1 | 0.46 | 0.50 | 0.00 | 0.00 | 0.00 | 1.00 | 1.00 | ▇▁▁▁▇ |
| female | 0 | 1 | 0.52 | 0.50 | 0.00 | 0.00 | 1.00 | 1.00 | 1.00 | ▇▁▁▁▇ |
| parental_schooling | 0 | 1 | 11.32 | 1.10 | 10.00 | 11.00 | 11.00 | 12.00 | 23.00 | ▇▁▁▁▁ |
| parental_lincome | 0 | 1 | 14.56 | 0.69 | 12.67 | 14.11 | 14.52 | 14.95 | 19.45 | ▂▇▁▁▁ |
| letter | 0 | 1 | 0.25 | 0.43 | 0.00 | 0.00 | 0.00 | 0.00 | 1.00 | ▇▁▁▁▂ |
| learnings | 0 | 1 | 10.06 | 1.28 | 5.39 | 9.18 | 10.04 | 10.96 | 14.62 | ▁▃▇▃▁ |
| year | 0 | 1 | 6.00 | 2.58 | 2.00 | 4.00 | 6.00 | 8.00 | 10.00 | ▇▇▃▇▇ |
| test_score | 0 | 1 | 2.37 | 0.72 | -0.57 | 1.88 | 2.33 | 2.84 | 5.05 | ▁▂▇▃▁ |
Hemos removido todos los missing.
3.3 Modificar los datos
Un último paso que haremos antes de hacer estadística decsriptiva es modificar los nombres de algunas columnas para que se vean bien en la tablas.
Vambos renombrar la variable preu a preuniversitario. Lo hacemos con rename() del paquete dplyr. Esta función tiene una sintaxix similar a filter()
# renombremos preu a preuniversitario
analysisdata<-rename(school_data_selected, summerschool=summercamp)
# use head to view the first 6 observations
head(analysisdata)# A tibble: 6 × 10
person_id school…¹ summe…² female paren…³ paren…⁴ letter learn…⁵ year test_…⁶
<dbl> <dbl> <int> <int> <int> <dbl> <dbl> <dbl> <int> <dbl>
1 1 5 0 1 10 12.9 0 10.2 6 1.64
2 1 5 0 1 10 12.9 0 10.2 2 1.09
3 1 5 0 1 10 12.9 0 10.2 3 1.91
4 1 5 0 1 10 12.9 0 10.2 4 2.07
5 1 5 0 1 10 12.9 0 10.2 7 2.38
6 1 5 0 1 10 12.9 0 10.2 8 2.03
# … with abbreviated variable names ¹school_id, ²summerschool,
# ³parental_schooling, ⁴parental_lincome, ⁵learnings, ⁶test_scoreEn un siguiente paso, vamos a transformar los puntajes en la prueba PSU y las notas a una variable que tenga media 0 y desviación estándar 1. Nosotros como chilenos cnocemos el sistema, pero un lector de la investigación probablemente no. Es por esto que es mejor trabajar con variables estandarizadas.
# Standardize test score
# Group analysisdata by year
analysisdata<-group_by(analysisdata,year)
# Create a new variable with mutate
analysisdata<-mutate(analysisdata, test_score=(test_score-mean(test_score))/sd(test_score))Esto tambien podemos hacerlo dentro de un pipe:
# Standardize test score
# Group analysisdata by year
analysisdata<-group_by(analysisdata,year)
# Create a new variable with mutate
analysisdata<-mutate(analysisdata, test_score=(test_score-mean(test_score))/sd(test_score))
# show mean of test_score
print(paste("Mean of test score:",mean(analysisdata$test_score)))[1] "Mean of test score: 1.66440088479111e-17"Podemos comprobar que efectivamente la media y desviación estándar corresponden a dichos valores:
# show mean of test_score
print(paste("Mean of test score:",mean(analysisdata$test_score)))[1] "Mean of test score: 1.66440088479111e-17"#show sd of test_score
print(paste("SD of test score:",sd(analysisdata$test_score)))[1] "SD of test score: 0.999872448979073"Notar que si aplicamos mean() o sd() o cualquier otra función matemática a columnas que tienen valores missing, también dará un valor missing a menos que usemos la opción na.rm=FALSE.
Ya estamos bien, ahora pasamos a conocer mejor nuestros datos con estadística descriptiva.
4. Estadística descriptiva
Hasta ahora, cargamos datos en diversos formatos (csv, dta y xlsx) los unimos, re-estructuramos el dataset, removimos valores missing y generamos algunas transformaciones. El siguiente paso es empezar a conocer nuestros datos. Para esto haremos tablas de estadísticas descriptivas y también algunos graficos descriptivos.
4.1 Tablas de estadística descriptiva
Hasta ahora, ya conocemos dos maneras de calcular estadísticas resumen:
sumary()de R base. Esta función en realidad funciona en muchos tipos de objetos de R y suele dar un bien resumen. Pero no en el formato de una tabla exportable a un documento latex, word o etc. que podamos presentar en nuestra investigación o resultados.skim()del paquete skimr
4.1.1 Tabla de estadísticas descriptivas “lista para llevar”
Una forma rápida de obtener una tabla de estadísticas descriptivas es con un primo de skim() del paquete modelsummary.
# load modelsummary
library("modelsummary")
# create a summary stat table in Latex format
analysisdata%>%
filter(year==2)%>%
select(female,starts_with("paren"),letter,summerschool,test_score)%>%
datasummary_skim( fmt="%.2f")| Unique (#) | Missing (%) | Mean | SD | Min | Median | Max | ||
|---|---|---|---|---|---|---|---|---|
| year | 1 | 0 | 2.00 | 0.00 | 2.00 | 2.00 | 2.00 | |
| female | 2 | 0 | 0.52 | 0.50 | 0.00 | 1.00 | 1.00 | |
| parental_schooling | 12 | 0 | 11.32 | 1.10 | 10.00 | 11.00 | 23.00 | |
| parental_lincome | 3486 | 0 | 14.56 | 0.69 | 12.67 | 14.52 | 19.45 | |
| letter | 2 | 0 | 0.25 | 0.43 | 0.00 | 0.00 | 1.00 | |
| summerschool | 2 | 0 | 0.46 | 0.50 | 0.00 | 0.00 | 1.00 | |
| test_score | 3486 | 0 | 0.00 | 1.00 | −4.13 | 0.02 | 3.71 |
Esta se puede exportar a varios formatos, como word o latex con el parámetro output=["ruta donde guardar la tabla"]. Primero hagámoslo en word:
# load modelsummary
library("modelsummary")
# create a summary stat table in Latex format
analysisdata%>%
filter(year==2)%>%
select(female,starts_with("paren"),letter,summerschool,test_score)%>%
datasummary_skim( fmt="%.2f",
histogram=FALSE, output="tab_summary_statistics.docx")[1] "tab_summary_statistics.docx"La guardó en la carpeta por default. Si quieremos que esté en nuestra carpeta de output, podemos usar el paquete here::
# load modelsummary
library("modelsummary")
# create a summary stat table in Latex format
analysisdata%>%
filter(year==2)%>%
select(female,starts_with("paren"),letter,summerschool,test_score)%>%
datasummary_skim( fmt="%.2f",
histogram=FALSE, output=here::here("output/tab_summary_statistics.docx") )[1] "/Users/melita/Documents/GitHub/ayudantiasccs/output/tab_summary_statistics.docx"Ahora la hacemos en formato latex:
# load modelsummary
#library("modelsummary")
# create a summary stat table in Latex format
#analysisdata%>%
# filter(year==2)%>%
# select(female,starts_with("paren"),letter,summerschool,test_score)%>%
# datasummary_skim( fmt="%.2f",
# histogram=FALSE, output=here::here("output/tab_summary_statistics.tex"))4.1.2 Tablas customizadas
Para customizar nuestra tabla aun más, podemos usar la función datasummary() tambien del pquete modelsummary. Esta función perimte que definamos una fórmula de la estructura de la tabla.
# creamos una tabla de estadísticas descriptivas resumen
datasummary(female+parental_schooling+
letter+test_score~Factor(summerschool)*(Mean+SD),
sparse_header = FALSE,
data=filter(analysisdata,year==2))| Mean | SD | Mean | SD | |
|---|---|---|---|---|
| female | 0.51 | 0.50 | 0.52 | 0.50 |
| parental_schooling | 11.03 | 0.74 | 11.67 | 1.34 |
| letter | 0.09 | 0.29 | 0.43 | 0.49 |
| test_score | −0.30 | 0.85 | 0.35 | 1.05 |
En este ejemplo: - Listamos las variables a incluir separadas por un + similar to in a female+parental_schooling+pa... - Usamos a ~ para separar la lista de variables en la f´romula - Usamos la formula Factor(summerschool)*(Mean+SD) para mostror la media y desviación estándar por separado para cada grupo creado por la variable preunivsersitario. - Usamos Factor() para indicarle a R que debería considerar preuniversitario como una variable binaria.
- También podríamos haber hecho esto al limpiar y procesar la base de datos. - Podemos tambien invertir el orden (Mean+SD)*Factor(summerschool), lo que entonces daría primero la media y devsicacion estándar y luego separar por los valores de preuniversitario.
- Usamos la opción
sparce_header=FALSEpara especificar que queremos incluir e preuniversitario como título.
4.1.3 Nombres de variables
Hasta ahora hemos utilizado los nombres de variables directamente en las tablas. Estó no es muy estético, podemos cambiarle el nombre directamente con espacios y mayúsculas en el nombre para darle un emjor aspecto. También es posible asignanrle una “label” o etiqueta cuando creamos la tabla, como lo vemos en el ejemplo:
# load modelsummary
library("modelsummary")
# create a summary stat table
datasummary((`Female`=female)+
(`Parental schooling (years)`=parental_schooling)+
(`Parental income (log)`=parental_lincome)+
(`Received reminder letter`=letter)+
(`Test Score`=test_score)~
(`Attended summer school`=Factor(summerschool))*
(Mean+SD),
sparse_header = FALSE,
data=filter(analysisdata,year==2))| Mean | SD | Mean | SD | |
|---|---|---|---|---|
| Female | 0.51 | 0.50 | 0.52 | 0.50 |
| Parental schooling (years) | 11.03 | 0.74 | 11.67 | 1.34 |
| Parental income (log) | 14.35 | 0.51 | 14.81 | 0.79 |
| Received reminder letter | 0.09 | 0.29 | 0.43 | 0.49 |
| Test Score | −0.30 | 0.85 | 0.35 | 1.05 |
4.1.4 Exportando nuestras tablas
Podemos exportar nuestras tablas a word o Latex podemos usar la expresión output="[nombre del archivo y ruta]". En datasummary() se ve así:
# load modelsummary
library("modelsummary")
# create a summary stat table
datasummary((`Female`=female)+
(`Parental schooling (years)`=parental_schooling)+
(`Parental income (log)`=parental_lincome)+
(`Received reminder letter`=letter)+
(`Test Score`=test_score)~
(Mean+SD+P25+P50+P75),
sparse_header = FALSE,
data=filter(analysisdata,year==2),
output = 'tab_descriptive_statistics.docx')4.2 Gráficos de estadística descriptiva
Vamos a usar principalmente la librería ggplot2 para crear nuestros gráficos.
4.2.1 Scatter plot (o gráfico de dispersión)
Nuestro primer grafico es un gráfico de dispersión. En este queremos ver como dos variables se relacionan en los datos. En estos podemos inlcuir curvas que describan la relación.
En este gráfico
- Iniciamos un objeto de grafico
ggplot()usando los datos analysisdata que ya procesamos y la vamos a filtrar solo para incluir el año 5. - Especificamos que parental_income sea el eje x y test_score el eje y en
aes() - Usamos
geom_point()para incluir los puntos que describen la dispersión. - Usamos
geom_smooth()para agregar una linea que describa la relación. - Usamos
theme()para darle formato a los elementos - Usamos
labs()para incorporar etiquetas a los ejes y al título.
# load ggplot2
library("ggplot2")
# create a scatter plot with a fitted line
ggplot(analysisdata%>%filter(year==5),
aes(x=parental_lincome,y=test_score))+
geom_point(alpha=0.1,size=0.85,color="#63a668")+
geom_smooth(color="#145c21") +
theme(panel.background = element_rect(fill="#ededed",color="#ededed"),
plot.background = element_rect(fill="#ededed",color="#ededed"),
panel.grid.major = element_line(colour="#a3a3a3",size=0.1))+
labs(x="Log(Parental Income)",y="Test Score (Mean=0,SD=1)", title="Test scores & Parental income")
4.2.2 Graficos de barras y boxplot
Del grafico anteriro podemos observar que los resultados de los test se correlacionan con ingreso parental. Esto no es una sorpresa. Veamos si tambien el asistir al preuniversitario se correlaciona con estas características individuales.
Primero, creamos un scatter plot de eduación de los padres y test score en el semestre antes del preuniversitario. Usemos el mismo código de arriba pero en lugar de ingreso, usamos educación parental. Segundo, creemos un gráfico de barras que muestre que el promedio de puntaje en la prueba diagnostico por asistencia al preuniversitario. Tercero, creamos un box plot de ingreso de los padres por asistencia al preuniversitario.
Nuestro codigo ahora tiene estos elementos.
- Creamos el objeto
ggplot()y cargamos la data y el tema. Este objeto es llamado un rawchart - Creamos 3 gráficos basados en rawchart, cada uno lo guardamos con un nombre.
- Para el de barras, usamos
geom_bar()y fijamosstat="summary", fun="mean"para decirle a R que cree un grafico de barras con la media de test_score
- usamos
labs()para decir los ejes y titulos. - Usamos
geom_boxplot()para crear el boxplot.
- Para el de barras, usamos
- Usamos el paquete patchwork para combinar varios gráficos en uno.
- usamos ggsave() para guardar el gráfico combinado como un archivo png
# Load patchwork
library("patchwork")
# Create raw chart element
rawchart<-ggplot(analysisdata%>%filter(year==4),x=as.factor(fill))+
theme_classic()
# Create bar chart of pre summer school test score and summer school
p1<-rawchart+
geom_smooth(aes(x=parental_schooling,y=test_score)) +
geom_point(aes(x=parental_schooling,y=test_score),alpha=0.1)+
labs(x="Parental schooling", y="Test Score Year 5")
# Create bar chart of pre summer school test score and summer school
p2<-rawchart+
geom_bar(aes(x=as.factor(summerschool),y=test_score),
stat="summary",fun="mean")+
labs(y="Test Score Year 5", x="Attended Summer School")
# Create bar chart of parental schooling and summer school attendance
p3<-rawchart+
geom_boxplot(aes(x=as.factor(summerschool),y=parental_lincome))+
labs(y="Parental Income (log)", x="Attended Summer School")
# Combine charts
p1/(p2+p3)
# Export chart
ggsave("fig1.png")Estos tres gráficos nos muestran que los puntajes en las pruebas se correlacionan con características de los padres (ambos scatters) y que aquellos que asistieron al preuniversitario tenían mejores puntaje sincluso antes de hacer el preuniversitario y que las caraterísticas de los padres se relaciona con la asistencia al preuniversitario.
Es decir, tenemos SESGO DE SELECCIÓN
4.1.3 Histogramas y gráficos de densidad
Comaperemos la distribución de los puntajes depués de asistir al preuniversitario por quienes fueron y los que no fueron a dicho reu. Crea,os un histograma con una linea que muestre las densidades estimadas.
# create a histogram and density chart
ggplot(filter(analysisdata,year==6),
aes(x=test_score,fill=as.factor(summerschool)))+
geom_histogram(aes(y=..density..),bins = 50,alpha=0.5,
position="identity",color="white")+
geom_density(alpha=0.0,size=1,show.legend= FALSE)+
theme_minimal()+
labs(y="Densidad",x="Puntaje en prueba de selección (estandarizado)",fill=" ")+
scale_fill_brewer(palette="Set2",labels=c("No asistió","Asistío a preuniversitario"))+
theme(legend.position="top")
Claramente hay diferencias en sus resultados. Pero la estimación directa va a confundir cuanto de esto proviene del sesgo de selección y cuanto es el efecto real de ir al preu.
4.1.4 Correalograma
Muchas veces queremos saber que tán correlacionadas estan las variables en una muestra. Podríamos simplementa calular una tabla de correlaciones, pero el paquete [PerformanceAnalytics] tiene una función muy conveniente: chart.Correlation() que nos presenta un gráico con las correlaciones de a pares, su significancia estadística, gráficos de dispersión y distribución.
#Cargamos el paquete
library("PerformanceAnalytics")
correlations_graph <- analysisdata %>%
dplyr::filter(year==6) %>%
ungroup() %>%
dplyr::select(female, parental_schooling, parental_lincome, summerschool, letter ) #variables para el correalograma
chart.Correlation(correlations_graph, histogram=TRUE, pch=19)
III. Estimación e identifiación
Tenemos nuestro problema base claro, ahora usemos algunas estrategias de identificaión para encontrar el efecto de asitir a un preuniversitario en esta muestra.
Para esto, partiremos del mejor caso: contamos con un experimento y luego iremos variando la información disponible y la estrategia que es factible emplear.
1. Regresión lineal ingenua
Nuestro objetivo es estimar la relación el puntaje en la prueba de admisión a la universidad y participar del preuniversitario intensivo. Para esto podríamos pensar en un modelo que directamente incluya a etsos dos elementos.
\[ Puntaje_i = \beta_0 + \beta_1 I(PREU=0)_i + u_i \] El modelo estimado sería:
ols <- lm ( test_score ~ summerschool, data= analysisdata )
ols
Call:
lm(formula = test_score ~ summerschool, data = analysisdata)
Coefficients:
(Intercept) summerschool
-0.4273 0.9309 summary(ols)
Call:
lm(formula = test_score ~ summerschool, data = analysisdata)
Residuals:
Min 1Q Median 3Q Max
-3.7025 -0.5296 0.0601 0.5993 3.2100
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.42732 0.00680 -62.84 <2e-16 ***
summerschool 0.93095 0.01004 92.75 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.8858 on 31361 degrees of freedom
Multiple R-squared: 0.2153, Adjusted R-squared: 0.2152
F-statistic: 8603 on 1 and 31361 DF, p-value: < 2.2e-16Podemos hacer una mejor tabla usando el paquete stargazer
library(stargazer)
stargazer::stargazer(ols, type = "text", out="tabla1.tex")
================================================
Dependent variable:
----------------------------
test_score
------------------------------------------------
summerschool 0.931***
(0.010)
Constant -0.427***
(0.007)
------------------------------------------------
Observations 31,363
R2 0.215
Adjusted R2 0.215
Residual Std. Error 0.886 (df = 31361)
F Statistic 8,603.086*** (df = 1; 31361)
================================================
Note: *p<0.1; **p<0.05; ***p<0.01Stargazer es util, pero esta siendo reemplazado por otros paquetes, particularmente modelsummary. Para una refrencia extra recomiendo revisar este link
# Podemos hacer un mapeo de nombres de variables
cm <- c( "test_score" = "Test Score",
"summerschool" = "Pre universitario" )library(modelsummary)
modelsummary::modelsummary(ols, stars= TRUE, fmt=3,
estimate = "{estimate} ({std.error}){stars}", statistic=NULL, gof_omit = "AIC|BIC|Lik",
coef_map = cm)| Model 1 | |
|---|---|
| Pre universitario | 0.931 (0.010)*** |
| Num.Obs. | 31363 |
| R2 | 0.215 |
| R2 Adj. | 0.215 |
| F | 8603.086 |
| RMSE | 0.89 |
También permite graficar:
modelsummary::modelplot(ols, coef_omit = "intercept")
No se ve muy bien para solo 1 modelo. Otra opcion es con el paquete parameters (https://easystats.github.io/parameters/)
library(parameters)
m <- parameters::model_parameters(ols)
mParameter | Coefficient | SE | 95% CI | t(31361) | p
--------------------------------------------------------------------------
(Intercept) | -0.43 | 6.80e-03 | [-0.44, -0.41] | -62.84 | < .001
summerschool | 0.93 | 0.01 | [ 0.91, 0.95] | 92.75 | < .001Este modelo, como ya lo revisamos de la estadística descriptiva, podemos pensar que sufre de sesgo de selección. Es posible que los estudiantes con padres más involucrados, con más educación y mayor ingreso
Si solo fueran el problema la eduación de los padres y su ingreso, podríamos controlar por estos. Incluyámoslos al modelo en “cascada”. Una manera de hacerlo es con varios modelos en una lista.
#modelos
naive_ols_models <- list(
" (1)" = lm( test_score ~ summerschool , data= analysisdata),
"(2)" = lm ( test_score ~ summerschool + parental_lincome, data= analysisdata ),
"(3)" = lm ( test_score ~ summerschool + parental_schooling, data= analysisdata ),
"(4)" = lm ( test_score ~ summerschool + parental_lincome + parental_schooling, data= analysisdata ) )
# Mapeo de nombre de variables
cm <- c( "test_score" = "Test Score",
"summerschool" = "Pre universitario",
"parental_lincome" = "Parental Income (log)",
"parental_schooling" = "Parental schooling" )Tabla con stargazer:
library(stargazer)
stargazer::stargazer(naive_ols_models, type = "text", out="tabla1.tex")
=========================================================================================================================================
Dependent variable:
---------------------------------------------------------------------------------------------------------------------
test_score
(1) (2) (3) (4)
-----------------------------------------------------------------------------------------------------------------------------------------
summerschool 0.931*** 0.551*** 0.685*** 0.554***
(0.010) (0.008) (0.009) (0.008)
parental_lincome 0.824*** 0.889***
(0.006) (0.010)
parental_schooling 0.381*** -0.051***
(0.004) (0.006)
Constant -0.427*** -12.255*** -4.626*** -12.622***
(0.007) (0.087) (0.047) (0.097)
-----------------------------------------------------------------------------------------------------------------------------------------
Observations 31,363 31,363 31,363 31,363
R2 0.215 0.506 0.377 0.507
Adjusted R2 0.215 0.506 0.377 0.507
Residual Std. Error 0.886 (df = 31361) 0.703 (df = 31360) 0.789 (df = 31360) 0.702 (df = 31359)
F Statistic 8,603.086*** (df = 1; 31361) 16,082.240*** (df = 2; 31360) 9,495.218*** (df = 2; 31360) 10,769.820*** (df = 3; 31359)
=========================================================================================================================================
Note: *p<0.1; **p<0.05; ***p<0.01Tabla con modelsummary:
modelsummary::modelsummary(naive_ols_models, stars= TRUE, fmt=3,
estimate = "{estimate} ({std.error}){stars}", statistic=NULL, gof_omit = "AIC|BIC|Lik",
coef_map = cm)| (1) | (2) | (3) | (4) | |
|---|---|---|---|---|
| Pre universitario | 0.931 (0.010)*** | 0.551 (0.008)*** | 0.685 (0.009)*** | 0.554 (0.008)*** |
| Parental Income (log) | 0.824 (0.006)*** | 0.889 (0.010)*** | ||
| Parental schooling | 0.381 (0.004)*** | −0.051 (0.006)*** | ||
| Num.Obs. | 31363 | 31363 | 31363 | 31363 |
| R2 | 0.215 | 0.506 | 0.377 | 0.507 |
| R2 Adj. | 0.215 | 0.506 | 0.377 | 0.507 |
| F | 8603.086 | 16082.240 | 9495.218 | 10769.815 |
| RMSE | 0.89 | 0.70 | 0.79 | 0.70 |
Sin embargo, si creemos que lo que también correlacona es una inobservable como la preocupación / involucramiento de los padres aun preservan el problema de la autoselección, todavía tenemos presencia de sesgo por variables omitidad.
¿Qué hacer?
Vamos a revisra varias alternativas: - Experimental - Datos observacionales: 1. Método de diferencias (en datos experimentales)
[2. Logit y probit](#logit-y-probit)
[3. Variables Instrumentales]
[4. Regresión Discontinua]
[5. Diferencias en Diferencias]
[6. Efectos fijos, aleatorios y modelos jerárquicos]
[7. Matching]2. Método de diferencias (en datos experimentales)
Un primer caso, es que antes de asistir al preuniversitario se selecciona al azar un grupo de estudiantes para enviarles una carta de invitación que recuerda sus características y detalla como participar.
¿Podría ser este un buen experimento?
- Asignación aleatoria
- Tratamiento: recibir la carta -> no es directamente el preuniversitario. Nos permite ver el efecto causal de recibir la carta, pero no de asistir al preuniversitario. Vamos a usar una técnica para enfrentar etso posteriormente.
Por ahora, supongamos que conocer el efecto de la carta en el puntaje es una pregunta lo suficientemente interesante.
\(testscore_i=\beta_0+\beta_1Letter_i+u_i\).
Este tipo de modelos podemos estimarlos con el paquete linear modelos lm() o también, con feols() del paquete fixtest
2.1 Estimación OLS
2.1.2 Usando lm()
#modelos
ols_letter_experiment <- list(
" (1)" = lm( test_score ~ letter , data= analysisdata) )
# Mapeo de nombre de variables
cm <- c( "test_score" = "Test Score",
"summerschool" = "Pre universitario",
"parental_lincome" = "Parental Income (log)",
"parental_schooling" = "Parental schooling" )Tabla con stargazer:
library(stargazer)
stargazer::stargazer(ols_letter_experiment, type = "text", out="tabla1.tex")
===============================================
Dependent variable:
---------------------------
test_score
-----------------------------------------------
letter 0.136***
(0.013)
Constant -0.034***
(0.006)
-----------------------------------------------
Observations 31,363
R2 0.003
Adjusted R2 0.003
Residual Std. Error 0.998 (df = 31361)
F Statistic 108.471*** (df = 1; 31361)
===============================================
Note: *p<0.1; **p<0.05; ***p<0.01Podemos ver que los que recibieron la carta, tienen 0.13 más puntaje que los que no, manteniendo lo demás constante.
Si le creemos a la asignción aleatoria, entonces este efecto podría ser interpetado como casual.
Para saber si la signación aletaoria efectivamente separó a dos grupos comparables de tratados y controles lo que debemos hacer es probar el balance en obbservables. Esto quiere decir que al menos en las variables observables, efectivamente ambos grupos no son estadísticamente diferentes antes del tratamiento.
2.1 Balance en observables
2.1 Comparación de a pares
Una primera opción es realizar una prueba t en las variables observables en ambos grupos, usamos el comando t.tes() . Por ejemplo podriamos comparar el ingreso de los padres con parental_income~summerschool para probar si la media de test_score es significantemente diferente en los dos grupos especificados por la asistencia al preuniversitario.
# Realizamos la prueba T
t.test(parental_lincome~summerschool,data=analysisdata)
Welch Two Sample t-test
data: parental_lincome by summerschool
t = -60.203, df = 24016, p-value < 2.2e-16
alternative hypothesis: true difference in means between group 0 and group 1 is not equal to 0
95 percent confidence interval:
-0.4765320 -0.4464809
sample estimates:
mean in group 0 mean in group 1
14.35102 14.81252 Podemos notar que no hay diferencias significativas entre los grupos.
Esto tendríamos que realizarlo en todas las variables y se suele presentar en una tabla de balance
2.3 Tabla de balance
Una tabla de balance presenta la diferencia de medias para grupos tratados y de control. Lo que deseamos es que no haya diferencias significativas antes del tratamiento en los grupos.
Podemos hacer una tabla que se vea bien usando datasummary_balance() del paquete modelsummary. Cargamos tambien el paquete estimatr porque le permite a datasummary_balance() calcular e inluir pruebas t de a pares directamente en la tabla.
# Load libraries
library(modelsummary)
library(estimatr)
# Filter and modify data
testdata<-filter(analysisdata,year==5)
testdata<-ungroup(testdata)
testdata<-mutate(testdata,Treated=ifelse(letter==1,"Letter","No Letter"))
testdata<-select(testdata,female,parental_schooling,parental_lincome,test_score,Treated)
testdata<-rename(testdata,`Female`=female,
`Parental schooling (years)`=parental_schooling,
`Parental income (log)`=parental_lincome,
`Test Score`=test_score)
# Table with balancing test
datasummary_balance(~Treated,
data = testdata,
title = "Balance of pre-treatment variables",
notes = "Notes: This is a brilliant table!",
fmt= '%.5f',
dinm_statistic = "p.value")| Mean | Std. Dev. | Mean | Std. Dev. | Diff. in Means | p | |
|---|---|---|---|---|---|---|
| Female | 0.52033 | 0.49988 | 0.51355 | 0.49991 | -0.00677 | 0.73032 |
| Parental schooling (years) | 11.31359 | 0.97129 | 11.32952 | 1.14573 | 0.01593 | 0.69028 |
| Parental income (log) | 14.56467 | 0.66201 | 14.56289 | 0.70338 | -0.00178 | 0.94629 |
| Test Score | -0.01950 | 0.99990 | 0.00641 | 1.00014 | 0.02591 | 0.50968 |
| Notes: This is a brilliant table! |
Los reslutados de la tabla sugieren que no hay diferencias significativas entre tratados y controles ANTES del tratamiento. ¡ Muy bien!
Y como antes, podríamos directaente agregar la opción outpu=tab_balancing.docx para exportarla en word (o latex).
2.4 Aproximación gráfica
Podemos ver el efecto causal propuesto en la siguiente gráfica. En los años antes del tratamiento son iguales y en los años posteriores, las distribuciones se empiezan a separar.
# Load ggridges
library("ggridges")
# create a ggridges chart
ggplot(analysisdata,aes(y=as.factor(year),x=test_score,fill=as.factor(letter) ))+
geom_density_ridges( alpha = .7, scale=1.5,color = "white", from = -2.5, to = 2.5)+
theme_minimal()+
theme_ridges(grid = FALSE)+
scale_y_discrete(expand = c(0, 0)) +
scale_x_continuous(expand = c(0, 0)) +
scale_fill_brewer(palette="Set1",labels=c("No letter","Letter"))+
labs(x="Test Score",y="Year",fill=" ",
title="Test score distribution by reminder letter status & year")+
theme(legend.position="top",aspect.ratio=4/3,plot.title = element_text(hjust = 0.5))
3. Efectos fijos, aleatorios y modelos jerárquicos
Muchas veces, los datos estan agrupados en diferentes categorías o dimensiones.
En nuestro caso los estudiantes están agrupados en colegios y es esperable que dentro de un colegio los resultados de dichos estudiantes estén correlacionados. Esto también ocurre al tener datos temporales o cualquier forma de agrupación lógica de la información.
Podemos ignorar esta estructura en los datos y simplemente considerarlos como un conjunto de unidades comparables, es decir tratarlo como un pooled cross-section. Si las diferencias entre grupos no son sistemáticas, no deberíamos observar diferencias, sin embargo esto puede estar haciendo que caigamos en problemas como la paradoja de Simpson o el Mathew effect
Para profundizar recomiendo: https://ds4ps.org/pe4ps-textbook/docs/p-040-fixed-effects.html o el capítulo
En nuestro caso, tenemos por ejemplo 30 colegios en los datos.
Otro punto es considerar que no solo el efecto está modulado por la organziación de agrupación de la data, sino también los errores. Esto lo llamamos clusterizar.
3.1.2 Estimación usando lm()
Si cada grupo está afectando los resultados, calculemos el efecto el un OLS solo del año 5, un pooled en los años hasta el 5 y en este ultimo, incluyamos un efecto fijo por colegio y luego por año.
# No siempre necesitamos hacer una lista de
OLS_year5_data <-filter(analysisdata,year==5)
OLS_year5 = lm( test_score ~ letter , data= OLS_year5_data )
# Para el pooled, usemos todos los años antes del 5
POLS_data <-filter(analysisdata,year<=5)
Pooled_OLS = lm ( test_score ~ letter , data= POLS_data )
# Fixed effect por colegio
FE = lm ( test_score ~ letter + factor(school_id), data= POLS_data )
# Usemos stargazer para la tabla esta vez
stargazer( OLS_year5, Pooled_OLS, FE,
type = "text",
dep.var.labels = ("Test score"),
column.labels = c("Cross-Sectional OLS (year 5)", "Pooled OLS", "School FE"),
covariate.labels = c("Intercept", "letter"),
keep =c("Intercept", "letter"),
omit.stat = "all",
digits = 2, intercept.bottom = FALSE )
===========================================================
Dependent variable:
-------------------------------------------------
Test score
Cross-Sectional OLS (year 5) Pooled OLS School FE
(1) (2) (3)
-----------------------------------------------------------
Intercept -0.03 -0.01 -0.02
(0.04) (0.02) (0.02)
===========================================================
===========================================================
Note: *p<0.1; **p<0.05; ***p<0.013.1.2 Estimación usando feols()
En el caso de que deseemos hacer un modelamiento más complejo, como opr ejemplo incluir efectos fijos (lo revisaremos al final de la sesión) o calcular los errores estándar en cluster podemos usar la función feols() del paquete fixest.
El paquete fixest permite la inclusión de múltiples efectos fijos y clusterización de errores es, nosotros usaremos la función feols()
Esta función espera una fórmula de especificación de la forma:
y~x1+x2+..|fixed effects|IV specification,cluster=..
Por ejemplo, en ese caso queremos clusterizar los errores a nivel de colegio:
# Load packages
library(fixest)
# Select data
regdata<-filter(analysisdata,year==6)
# Regression
m1<-feols(test_score~letter+parental_lincome+female+parental_schooling, cluster="school_id",data=regdata)
# Summary of regression
summary(m1)OLS estimation, Dep. Var.: test_score
Observations: 3,481
Standard-errors: Clustered (school_id)
Estimate Std. Error t value Pr(>|t|)
(Intercept) -14.571682 0.331506 -43.955950 < 2.2e-16 ***
letter 0.206422 0.026758 7.714325 1.6593e-08 ***
parental_lincome 1.019714 0.035335 28.858689 < 2.2e-16 ***
female -0.014779 0.020544 -0.719393 4.7766e-01
parental_schooling -0.028398 0.022283 -1.274430 2.1263e-01
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
RMSE: 0.725145 Adj. R2: 0.473408Quisieramos compararlo con el modelo sin dicha clusterización, podriamos hacer una mejor tabla. Aprovechamos la función modelsummary() y:
Corremos ambas regresiones y guardamos los modelos en una lista. La primera es sin ajustar los errores, la segunda con cluster, la tercera con un efecto fijo por colegio.
Usamos
modelsummary()para el output de las regresiones.Especificamos que entre paréntesis esté el error estándar y que omita la constante en la tabla.
Definimos el formato de los números con
fmt=...
# Load flextable and modelsummary
library(modelsummary)
# Select data
regdata<-filter(analysisdata,year==6)
# Regressions
models<-list(
m1<-lm(test_score~letter+parental_lincome+female+parental_schooling ,data=regdata),
m2<-feols(test_score~letter+parental_lincome+female+parental_schooling, cluster="school_id",data=regdata),
m3<-feols(test_score~letter+parental_lincome+female+parental_schooling|school_id, cluster="school_id",data=regdata)
)
# Generate table
modelsummary(models, stars = TRUE, statistic = 'std.error',
fmt= '%.4f',
coef_omit= '(Intercept)')| Model 1 | Model 2 | Model 3 | |
|---|---|---|---|
| letter | 0.2064*** | 0.2064*** | 0.2110*** |
| (0.0285) | (0.0268) | (0.0276) | |
| parental_lincome | 1.0197*** | 1.0197*** | 1.0222*** |
| (0.0299) | (0.0353) | (0.0363) | |
| female | −0.0148 | −0.0148 | −0.0100 |
| (0.0246) | (0.0205) | (0.0200) | |
| parental_schooling | −0.0284 | −0.0284 | −0.0310 |
| (0.0187) | (0.0223) | (0.0227) | |
| Num.Obs. | 3481 | 3481 | 3481 |
| R2 | 0.474 | 0.474 | 0.480 |
| R2 Adj. | 0.473 | 0.473 | 0.475 |
| R2 Within | 0.473 | ||
| R2 Pseudo | |||
| AIC | 7653.2 | 7651.2 | 7669.0 |
| BIC | 7690.1 | 7682.0 | 7878.3 |
| Log.Lik. | −3820.589 | −3820.589 | −3800.520 |
| F | 783.133 | ||
| RMSE | 0.73 | ||
| Std.Errors | by: school_id | by: school_id | |
| FE: school_id | X | ||
| + p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001 |
Podemos obsercar que no cambian los estimadores entre el modelo 1 y 2, pero si el error. Podemos ver que en el modelo 3 también cambia el estimador.
3.1.3 Elementos clave de los efectos fijos:
Los elementos centrales de efectos fijos son:
Usar efectos fijos nos permite controlar por todas las características invariantes en la unidad de analisis.
Asumimos que algunas de las características de la entidad que estamos estudiando afectan nuestra prediccion y necesitamos controlar por estas. Notar que efectos fijos nos permite entonces, capturar todas las caracteristicas que no varian en el tiempo en dicha entidad reduciendo indirectamente el sesgo de variables omitidas.
Importante saber que no podemos estudiar el efecto de estas caracter+pisticas que no varian en el tiempo usando efectos fijos. Por ejemplo, no podríamos estudiar ni el impacto del sexo, ni de la educación de los padres en este ejemplo. Como dichas variables no cambian en el tiempo, la constante que agregamos en el efecto fijo absorve dicho efecto tambien, pues son perefectaente colineares.
Por esto, es importante tener en cuenta que los efectos fijos estan diseñados para estudiar características que cambian en el tiempo para una entidad.
_ finalmente. efectos fijos tambien asume que cada entidad es única, tal que el error de dicha entidad es constante y no se correlaciona con el error de otra entidad. Es decir, asume independencia entre entidades.
3.2 Random vs fixed effects
CUn supuesto subyacente de la implementación de efectos fijos es que no hay correlación entre los errores de las entidades de agrucación y las constantes, es decir, que cada colegio es indepediente del otro (en nuestro caso).
Este supuesto se puede testear estadísticamente mediante un modelo aletrnativo que, si bien calcula agrupando los datos permite que los efectos sean aleatorios, es decir calculamos asumiendo que las diferencias entre las entidades afectan a la variable dependiente y que hau correlación entre el error de las entidades y sus constantes. Podemos usar la función plm() del paquete plm especificando model=“random”.
Necesitamos que el panel esté bien conformado para poder aplicar este paquete, es decir tener una observación por unidad dimensionada y variables que cambien intra y entre grupos.
#regdata <- analysisdata %>%
# group_by(school_id, year)%>%
# dplyr::summarize( test_score = mean(test_score), summerschool= mean(letter), female= mean(female), suma=n() ) %>%
# ungroup()
#head(regdata)# cargamos la librería
#library(plm)
# Para especificar un modelo correctaente, necesitamos variables entre e intra entidades. EN este caso, agruparemos la data por colegio y año.
#regdata <- analysisdata %>%
# group_by(school_id, year) %>%
# summary( test_score = mean(test_score), summerschool=test_score ) %>%
# ungroup()
#Fixed effects
#fixed <- plm( test_score ~ summerschool, data=regdata, index=c("school_id", "year"), model="within" )
#Random effects
#random <- plm( test_score ~ summerschool, data=regdata, index=c("school_id", "year"), model="random" )Podemos comparar ambos modelos mediante la prueba de Hausman para decidir que modelo es mejor. Generalmente se prefiere efecto fijo ya que es más eficiente, si ambos están correctos.
Algunos otros elementos a tener en cuenta: - Los efectos aleatorios permiten icnluir elementos que son invariables en el tiempo dado que asume que el error de la entidad no está correlacionado con el predictor. - Los efectos aleatorios están influenciados por sesgo de variable omitida, ya que no estamos incluyendo directamente las dummies.
3. Logit y probit
Un segundo elemento es revisar que ocurre con probabilidad de entrar al preuniversitario dado que se recibió la carta.
En este caso, el modelo a estimar tendría como outcome de interés a una variable dicotómica.
\[ preuniversitario_i = g(letter_i)+u_i \]
3.1 Probabilidad lineal
Una primera aproximación sería sencillamente estimar un modelo lineal:
\[ preuniversitario_i = \beta_0+ \beta_1 letter_i + \gamma' X_i+u_i \] Ya sabemos como hacerlo, simplemente usamos lm() o felm(). Aprovechemos esta oportunidad para agregar una fila con la media de la variable dependiente, este es genralemnte un elemento importante a incluir cuando estamos rabajando con variables dicotómicas ya que pone en contexto el análisis. Esto lo hacemos con el parámetro add_rows=...
# Estimate LPM
models<-list(
m1<-feols(summerschool ~letter, cluster="school_id",data=regdata),
m2<-feols(summerschool ~letter+parental_schooling+parental_lincome+female,cluster="school_id",data=regdata)
)
# Store the mean of dependent variable in a data frame
added_stats<-tibble("Mean of Dep. ",m1=mean(regdata$summerschool),m2=mean(regdata$summerschool))
# Generate table
modelsummary(models, stars = TRUE,statistic = 'std.error',
fmt= '%.4f',add_rows = added_stats,
coef_omit= '(Intercept)', output = 'flextable')
| Model 1 | Model 2 |
letter | 0.4451*** | 0.4457*** |
(0.0134) | (0.0134) | |
parental_schooling | 0.0353** | |
(0.0120) | ||
parental_lincome | 0.1925*** | |
(0.0228) | ||
female | 0.0086 | |
(0.0104) | ||
Num.Obs. | 3481 | 3481 |
R2 | 0.148 | 0.260 |
R2 Adj. | 0.148 | 0.260 |
R2 Within | ||
R2 Pseudo | ||
AIC | 4474.4 | 3989.4 |
BIC | 4486.7 | 4020.2 |
Log.Lik. | -2235.192 | -1989.700 |
Std.Errors | by: school_id | by: school_id |
Mean of Dep. | 0.4591 | 0.4591 |
+ p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001 | ||
Podemos ver que recubir la carta aumenta la probabilidad de asistir al preuniversitario en 44%, lo que es un incremento de casi un 100% considerando que la media de dicha variable era 0.46.
A pesar de su gran interpretabilidad, los modelos de probabilidad lineal tienen varios problemas: predicciones pueden salir del rango [0,1] y que son heterocedasticos por definición.
Una alternativa a este modelamiento es usar una forma funciónal que asegure que la predicciones estén en el rango [0,1] es con logit o probit.
3.2 Logit y probit
Para estimar logit usamos glm() para estimar modelos logit y probit. Espeificamos la función en family=binomial(link="probit") o family=binomial(link="logit").
Estimemos con logit:
# Estimate a binary outcomes model using a logit
logit_results <- glm(summerschool ~letter, data = regdata, family = binomial(link="logit"))
# Print the results
summary(logit_results)
Call:
glm(formula = summerschool ~ letter, family = binomial(link = "logit"),
data = regdata)
Deviance Residuals:
Min 1Q Median 3Q Max
-1.7781 -0.9267 -0.9267 1.4508 1.4508
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.62298 0.04098 -15.20 <2e-16 ***
letter 1.97333 0.09377 21.04 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 4802.3 on 3480 degrees of freedom
Residual deviance: 4265.4 on 3479 degrees of freedom
AIC: 4269.4
Number of Fisher Scoring iterations: 4El problema es que estos parámetros no son directamente interpretables porque son las pendientes en la función logistica. Para esto, calculamos los efectos marginales.
3.3 Efectos marginales
Usamos la librería margins para obtener los efectos marginales promedio. Simplemente usamos la función margins() en un objeto que tenga la estimación por probit o logit.
# Load margins package
library("margins")
# Select data
regdata<-filter(analysisdata,year==6)
## Estimate a binary outcomes model using a logit and probit
logit_results <- glm(summerschool ~letter+parental_lincome+female+parental_schooling, data = regdata, family = binomial(link="logit"))
probit_results <- glm(summerschool ~letter+parental_lincome+female+parental_schooling, data = regdata, family = binomial(link="probit"))
# Estimate linear probability model
lpm<-feols(summerschool ~letter+parental_schooling+parental_lincome+female, cluster="school_id",data=regdata)
# Compute marginal effects
mfx_logit<- margins(logit_results)
mfx_probit<- margins(probit_results)
# Include in table
modelsummary(list("LPM"=lpm,"Logit"=mfx_logit,"Probit"=mfx_probit))| LPM | Logit | Probit | |
|---|---|---|---|
| (Intercept) | −2.859 | ||
| (0.237) | |||
| letter | 0.446 | 0.408 | 0.414 |
| (0.013) | (0.013) | (0.013) | |
| parental_schooling | 0.035 | 0.046 | 0.048 |
| (0.012) | (0.012) | (0.012) | |
| parental_lincome | 0.193 | 0.189 | 0.183 |
| (0.023) | (0.017) | (0.017) | |
| female | 0.009 | 0.009 | 0.009 |
| (0.010) | (0.014) | (0.014) | |
| Num.Obs. | 3481 | 3481 | 3481 |
| R2 | 0.260 | ||
| R2 Adj. | 0.260 | ||
| R2 Within | |||
| R2 Pseudo | |||
| AIC | 3989.4 | 3778.4 | 3785.3 |
| BIC | 4020.2 | 3809.2 | 3816.0 |
| Log.Lik. | −1989.700 | −1884.212 | −1887.635 |
| Std.Errors | by: school_id |
Acá podemos ver una tabla con los tres modelos y compararlos.
4. Variables Instrumentales
Volviendo a nuestra pregunta de investigación inicial queríamos conocer elefecto de asistir al preuniversitario en los puntajes, no de la carta.
Sin embargo, sabemos que la carta se asignó aleatoriamente y que afecta a la probabilidad de asitir al preuniversitario. En otras palabras: podemos usar la carta como instrumento
Es más, podríamos no realizar ninguna estimación y recurrir a lo ya obtenido:
- Sabemos que el estimador de la forma reducida: \[ testscore_i = \beta_0 + \beta_1 letter_i + u_i \] es 0.2 SD.}
- En el punto anterior calculamos la primera etapa: \[ preuniversitario_o = \alpha_0 + \alpha_1 letter_i + \mu_i\] que corresponde aproximadamente a 0.4.
Entonces, el Estimador de Wald de variables instrumentales es: IV=0.2/0.4 =0.5
Podríamos obtener este estimador directamente ahora desde los datos. Vamos a revisar dos métodos.
4.2 IV usando ivreg()
Primero usemos la función del paquete AER (Applied Economics with R). La sintaxis es: - ivreg(outcome equation| exogenous variables, data)
# Load the AER package
library(AER)
# Select data
regdata<-filter(analysisdata,year==6)
# Estimate IV specification with
summary(ivreg(test_score~summerschool+female+parental_lincome+parental_schooling|
female+parental_lincome+parental_schooling+letter,data=regdata))
Call:
ivreg(formula = test_score ~ summerschool + female + parental_lincome +
parental_schooling | female + parental_lincome + parental_schooling +
letter, data = regdata)
Residuals:
Min 1Q Median 3Q Max
-2.334605 -0.443619 0.007695 0.435387 2.560634
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -13.24728 0.30615 -43.271 < 2e-16 ***
summerschool 0.46319 0.05723 8.093 7.96e-16 ***
female -0.01878 0.02203 -0.853 0.394
parental_lincome 0.93053 0.02898 32.107 < 2e-16 ***
parental_schooling -0.04476 0.01687 -2.654 0.008 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.649 on 3476 degrees of freedom
Multiple R-Squared: 0.5793, Adjusted R-squared: 0.5788
Wald test: 979.1 on 4 and 3476 DF, p-value: < 2.2e-16 Observamos que el estimador es cercano a 0.5.
4.2 IV usando feols()
Podemos usar nuevamente feols() e incluir variables instrumentales. Además felm() también permite estimar especificaciones. En este caso no incluirmos efecto fijo, así que le poneos 0.
# Estimate IV specification with felm
m1<-feols(test_score~parental_lincome+female+parental_schooling| # Outcome eq.
0| # Fixed effects
summerschool~letter # First stage
,cluster="school_id" # Cluster var
,data=regdata)
# Summary of results
summary(m1)TSLS estimation, Dep. Var.: test_score, Endo.: summerschool, Instr.: letter
Second stage: Dep. Var.: test_score
Observations: 3,481
Standard-errors: Clustered (school_id)
Estimate Std. Error t value Pr(>|t|)
(Intercept) -13.247279 0.326357 -40.591334 < 2.2e-16 ***
fit_summerschool 0.463188 0.058299 7.945045 9.2096e-09 ***
parental_lincome 0.930534 0.030648 30.362006 < 2.2e-16 ***
female -0.018783 0.018942 -0.991616 3.2959e-01
parental_schooling -0.044757 0.019024 -2.352669 2.5636e-02 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
RMSE: 0.648514 Adj. R2: 0.578825
F-test (1st stage), summerschool: stat = 699.1, p < 2.2e-16 , on 1 and 3,476 DoF.
Wu-Hausman: stat = 40.7, p = 1.96e-10, on 1 and 3,475 DoF.Los coeficientes son iguales que antes, pero los errores estándar cambiaron ya que los clusterizamos al nivel del colegio. Notemos que feols() nos indica que el F-stat del isntrumento excluido es 687 (cumpliendo la condición de relevancia del isntrumento).
4.3 Tabla de reporte IV
Terminemos la sección de IV con una tabla que muestre la correlación pira entre la asistencia al preuniversitario y puntajes (ingenuo), la forma reducida, la primera etapa y la estimación de variable instrumental
Notemos que usaremos coef_map= c( ) para renombrar las variables en la tabla, espeficicar los coeficientes a incluiry permitir que dos coeficientes con diferente nombre aparezcan en la misma fila (preuniversitario y fit_preuniversitario)
# Estimate OLS
OLS<-feols(test_score~summerschool+parental_lincome+female+parental_schooling,cluster="school_id",data=regdata)
# Estimate reduced form
RF<-feols(test_score~letter+parental_lincome+female+parental_schooling,cluster="school_id",data=regdata)
# Estimate first stage
FS<-feols(summerschool~letter+parental_lincome+female+parental_schooling,cluster="school_id",data=regdata)
# Estimate IV specification
IV<-feols(test_score~parental_lincome+female+parental_schooling| # Outcome eq.
0| # Fixed effects
summerschool~letter # First stage
,cluster="school_id" # Cluster var
,data=regdata)
# Combine results
IVresults<-list("OLS"=OLS,"RF"=RF,"FS"=FS,"IV"=IV)
# Coefficients
cm <- c('fit_summerschool' = 'Summer School',
'summerschool' = 'Summer School','parental_lincome' = 'Parental Income',
'letter' = 'Reminder letter', "female"="Female","parental_schooling"="Parental Schooling")
# Output Table
modelsummary(IVresults, stars = TRUE,statistic = 'std.error',
fmt= '%.4f',coef_map=cm, output = )| OLS | RF | FS | IV | |
|---|---|---|---|---|
| Summer School | 0.7855*** | 0.4632*** | ||
| (0.0160) | (0.0583) | |||
| Parental Income | 0.8671*** | 1.0197*** | 0.1925*** | 0.9305*** |
| (0.0278) | (0.0353) | (0.0228) | (0.0306) | |
| Reminder letter | 0.2064*** | 0.4457*** | ||
| (0.0268) | (0.0134) | |||
| Female | −0.0224 | −0.0148 | 0.0086 | −0.0188 |
| (0.0185) | (0.0205) | (0.0104) | (0.0189) | |
| Parental Schooling | −0.0551** | −0.0284 | 0.0353** | −0.0448* |
| (0.0164) | (0.0223) | (0.0120) | (0.0190) | |
| Num.Obs. | 3481 | 3481 | 3481 | 3481 |
| R2 | 0.602 | 0.474 | 0.260 | 0.579 |
| R2 Adj. | 0.602 | 0.473 | 0.260 | 0.579 |
| R2 Within | ||||
| R2 Pseudo | ||||
| AIC | 6678.6 | 7651.2 | 3989.4 | 6873.6 |
| BIC | 6709.3 | 7682.0 | 4020.2 | 6904.4 |
| Log.Lik. | −3334.286 | −3820.589 | −1989.700 | −3431.801 |
| Std.Errors | by: school_id | by: school_id | by: school_id | by: school_id |
| + p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001 |
La primera columna nos muestra que si ignoramos el problema de selección y etsimamos la regresión directamente entre el puntaje y la asistencia al preuniversitario (condicional en los controles) nos dá un efecto promedio de 0.78 SD. Como sería esperable, este efecto está sobrestimado pueso incluye el sesgo por selección, lo que podemos comprobar al comprarlo con el estimador de IV de 0.47SD.
Las estadísticas descriptivas que estudiamos, sugieren que quienes asistieron al preuniversitario ya tenían mejores puntajes antes del preuniversitario y que sus padres tenían mayores ingresos y años de educación. Entonces, el coeficiente en la primera columna no solo captura el efecto rpopio del preuniversitario, sino todas estas caracteristicas inobservables correlacionadas positivamente con los puntajes y con la asistencia al preuniversitario.
6. Diferencias en Diferencias
La estrategia que tomamos con IV es bastante convincente, porque se basa en un exterimento aleatorio (RCT) y el coeficiente “ingenuo” era mayor que el de IV, como era esperable.
¿Pero que hacemos si no contamos con una carta? Imaginemos que no tenemos la carta, pero si contamos con mediciones antes y despues del preuniversitario (como lo hacemos). Podemos usar la estrategia de DiD .
Estaría dada por la siguiente especificación:
\[test_score=\beta1+\beta_2summerschool_i+\beta_3after_i+\beta_4after_i\times summerschool+\gamma'X_i+u_i\]
en el cual after es una variable dictotómica que considera el periodo después de asistir al preuniversitario.
6.1 Tendencias paralelas
Un supuesto clave en la estrategia de diferencias en diferencias es que, en ausencia del tratamiento los grupos habrian estado evolucionando “en paralelo”. Con esto argumentamos que la evolución en los controles son un buen contrafactual a la evolución de los tratados.
# Collapse data on the year X summer school level
df<-group_by(analysisdata,year,summerschool)
# Calculate the mean test score by year
df<-summarise(df,test_score=mean(test_score))
# Create a chart
ggplot(df,aes(x=year,y=test_score,color=as.factor(summerschool)))+
geom_line(size=1)+
theme_classic()+
theme(legend.position = "top")
Observamos que previo al preuniversitario iban “mas o menos” en paralelo. Pero después del tratamiento, estos caen para lso controles. Eto e sun poco preocupante, pero fuera de eso, se ven similares.
Generalmente, no contamos con tantas mediciones para tener una evaluación empirica de las tendencias.
6.2 Estimación modelo de diferencias en diferencias.
Primero, en un caso más similar, usemos solo las estimaciones justo antes y una despues del preuniversitario usando la función lm().
Primero, creamos la variable dummy de después usando mutate() y ifelse().
# Select year 5 and year 6 observations
df<-analysisdata%>%filter(year%in%c(5,6))
# Generate an after indicator and the interaction between after and summer school
df<-mutate(df,after=ifelse(year==6,1,0),afterXsumerschool=after*summerschool)
# Estimate the DiD regression and save the results in an element called m1
m<-lm(test_score~summerschool+after+afterXsumerschool,data=df)
# Print the results
summary(m)
Call:
lm(formula = test_score ~ summerschool + after + afterXsumerschool,
data = df)
Residuals:
Min 1Q Median 3Q Max
-3.5022 -0.5325 0.0603 0.6072 3.3260
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.18933 0.02081 -9.098 <2e-16 ***
summerschool 0.41230 0.03071 13.426 <2e-16 ***
after -0.33844 0.02942 -11.502 <2e-16 ***
afterXsumerschool 0.73736 0.04343 16.980 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.9027 on 6957 degrees of freedom
Multiple R-squared: 0.1853, Adjusted R-squared: 0.1849
F-statistic: 527.4 on 3 and 6957 DF, p-value: < 2.2e-16Bien, observamos que el estiamodr de DIF es de 0.75SD. Observamos que es mayor que el que obtuvimos en IV y menor que el de OLS, pero por poco.
Veamos que ocurre si ahora usamos todos los datos disponibles. Podemos usar la función feols() para estimar con todos los años (también podriamos usar lm()). Este es el enfoque tradiciona, con una variable que indica después del tratamiento y la interacción entre después y tratamiento y una dummy por tratamiento. En el segundo modelo ponemos la otra opción, que es agregar un efecto fijo temporal en bez de la dummy de después. En el tercer modelo, agregamos controles.
# Modify data
df<-mutate(analysisdata,after=ifelse(year>5,1,0),afterXsummerschool=after*summerschool)
df<-rename(df,school=summerschool)
# Estimate DiD without controls and fixed effects
m1<-feols(test_score~school+afterXsummerschool+after,cluster="school_id",data=df)
# Estimate DiD without controls
m2<-feols(test_score~school+afterXsummerschool|school_id+year,cluster="school_id",data=df)
# Estimate DiD with controls
m3<-feols(test_score~afterXsummerschool+school+parental_lincome+female+parental_schooling
|school_id+year,cluster="school_id",data=df)
# Leave year 5 out to avoid dip
m4<-feols(test_score~afterXsummerschool+school+parental_lincome+female+parental_schooling
|school_id+year,cluster="school_id",data=filter(df,year!=5))
# Coefficients
cm <- c('afterXsummerschool' = 'Summer School X After',
'school' = 'Summer School',
'after'='After')
# Output Table
modelsummary(list(m1,m2,m3,m4), stars = TRUE,statistic = 'std.error',
statistic_override = vcov,
fmt= '%.4f',coef_map=cm, output = 'gt'
)| Model 1 | Model 2 | Model 3 | Model 4 | |
|---|---|---|---|---|
| Summer School X After | 0.6010*** | 0.6011*** | 0.6010*** | 0.5395*** |
| (0.0125) | (0.0125) | (0.0125) | (0.0132) | |
| Summer School | 0.5970*** | 0.5924*** | 0.2153*** | 0.2836*** |
| (0.0345) | (0.0344) | (0.0208) | (0.0211) | |
| After | -0.2758*** | |||
| (0.0077) | ||||
| Num.Obs. | 31363 | 31363 | 31363 | 27883 |
| R2 | 0.237 | 0.244 | 0.533 | 0.549 |
| R2 Adj. | 0.237 | 0.243 | 0.533 | 0.548 |
| R2 Within | 0.236 | 0.528 | 0.544 | |
| R2 Pseudo | ||||
| AIC | 80502.3 | 80296.3 | 65185.0 | 56995.6 |
| BIC | 80535.7 | 80630.5 | 65544.2 | 57341.5 |
| Log.Lik. | -40247.160 | -40108.165 | -32549.481 | -28455.818 |
| FE: school_id | X | X | X | |
| FE: year | X | X | X | |
| + p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001 | ||||
Podemos ver que el efecto de tratamiento es cercano a 0.6SD, lo que se acerca más al de IV y se aleja aun más que del OLS ingenuo.
La última columna (4) excluye el año 5 donde veiamos la desviación en las tendencias. Esto nos da un estimador aun más cerca del de IV.
8. Matching
Otra alternativa de análisis es mediante el matching. EN este caso, buscaremos para cada individuo que fue al preuniversitario algun compañero similar en observables que no haya asisitido.
Hay dos formas de hacerlo, buscar alguien que sea exactamente igual (exact matching) o que la probabilidad de asistir al preuniversitario basado en un set de observables sea igual (propensity score matching). Vamos a ver ambos.
Una revisión completa de matching está acá: ds4ps
La estrategia de matching sigue los siguientes dos pasos:
1 - Para cada tratado encontramos a su “clon”. Para esto hay dos alternativos, usamos uno (o varios) match exactos o usamos el propensity score matching. 2 - Usando los “clones” calculamos calculamos el ATE.
Lo primero que deberíamos hacer entonces es obtener los clones. Vamos a trabajar con el paquete MatchIT (Muy recomendable revisar su documentación)
#cargar el paquete
library("MatchIt")8.1 Exact matching
Un exact match, lo que hace es identificar uno ( varios) contrafactuales que se encuentren más cercanos a cada individuo tratado en la muestra.
Lo que hacemos es definir en que variables vamos a considerar el match. En este caso lo haremos por colegio, sexo, educacion de los padres e ingreso de los padres.
# Definamos un matching con todo por defecto.
a <- MatchIt::matchit( summerschool ~ school_id + female + parental_schooling + parental_lincome, data=analysisdata )
aA matchit object
- method: 1:1 nearest neighbor matching without replacement
- distance: Propensity score
- estimated with logistic regression
- number of obs.: 31363 (original), 28792 (matched)
- target estimand: ATT
- covariates: school_id, female, parental_schooling, parental_lincomeObservemos que lo que obtenemos es un objeto de matchit usando por edefto el vecino más cercano sin reemplazo.
summary(a)
Call:
MatchIt::matchit(formula = summerschool ~ school_id + female +
parental_schooling + parental_lincome, data = analysisdata)
Summary of Balance for All Data:
Means Treated Means Control Std. Mean Diff. Var. Ratio
distance 0.5226 0.4051 0.6197 2.1851
school_id 15.5788 15.7210 -0.0161 1.0578
female 0.5251 0.5074 0.0354 .
parental_schooling 11.6746 11.0281 0.4830 3.2933
parental_lincome 14.8125 14.3510 0.5854 2.3482
eCDF Mean eCDF Max
distance 0.1864 0.3147
school_id 0.0110 0.0331
female 0.0177 0.0177
parental_schooling 0.0536 0.2586
parental_lincome 0.1849 0.3127
Summary of Balance for Matched Data:
Means Treated Means Control Std. Mean Diff. Var. Ratio
distance 0.5226 0.4299 0.4886 2.4397
school_id 15.5788 15.7006 -0.0138 1.0616
female 0.5251 0.5169 0.0163 .
parental_schooling 11.6746 11.1325 0.4050 3.4192
parental_lincome 14.8125 14.4520 0.4574 2.7339
eCDF Mean eCDF Max Std. Pair Dist.
distance 0.1358 0.2825 0.4889
school_id 0.0116 0.0337 1.1368
female 0.0081 0.0081 1.0056
parental_schooling 0.0449 0.2281 0.4797
parental_lincome 0.1365 0.2806 0.4948
Sample Sizes:
Control Treated
All 16967 14396
Matched 14396 14396
Unmatched 2571 0
Discarded 0 0Podemos graficar el balance
plot(summary(a))
8.2 Propensity score matching
From the package MatchIt we will use the function matchit to create the propensity scores and the new dataset with the matched observation.
Note the function code. First, we specificy the logit model to calculate the propensity scores as we did before in section @ref:creating-matching-score. We then specify the method to match the scores. We will use the nearest neighbor method (see @ref:nearest-neighbor-matching-and-caliper-distance). The default method is “withour replacement”.
#m1 = matchit(df$school ~ df$income + df$occ_score + df$y_educ, method = "nearest", data = df)We can see the outcome of the process by using summary.
Have a look at the column “Mean Diff” and compare results in the “Summary of balance for all data” vs “Summary of balance for matched data”. You can see that the mean difference is drastically reduced in the matched data.
We can also look at the “Sample sizes” summary at the very end of the output. We started with 8841 control cases and 1159 treated cases. The algorithm matched all treated cases with only one control case, discarding 7682 observations.
#summary(m1)Once the algotrithm has calculated and matched the propensity scores, we can create a new dataset that contains all matched observations using the function match.data.
#m1data = match.data(m1)ow that we have the dataset with the matched observations, we can have another look at the covariates and compare the treated and untreated cases.
Para esto, hacemos una tabla de balance como lo hicismos anteriormente:
#t.test(m1data$y_educ[m1data$school==1], m1data$y_educ[m1data$school==0], paired = T)We can see this graphically by looking at the dentisy plots for each independent variables in the treated and untreated group.
#m1data$school = as.factor(m1data$school2)
#
# ggplot(m1data, aes(x = m1data$y_educ, fill = m1data$school2)) +
# geom_density(position="identity", alpha=0.6) +
# scale_fill_brewer(palette="Dark2") +
# ggtitle("Density plot of years of education, by group") +
# scale_x_continuous(name = "Years of education",
# breaks = seq(14, 25, 2),
# limits=c(14, 25)) +
# scale_y_continuous(name = "Density") +
# theme(plot.title = element_text(size = 14, face = "bold"),
# text = element_text(size = 12)) +
# guides(fill=guide_legend(title=NULL))+
# scale_fill_discrete(labels=c("Control - Public School", "Treatment - Private School"))# ggplot(m1data, aes(x = m1data$income, fill = m1data$school2)) +
# geom_density(position="identity", alpha=0.6) +
# scale_fill_brewer(palette="Dark2") +
# ggtitle("Density plot of income, by group") +
# scale_x_continuous(name = "Income",
# breaks = seq(11, 13, 0.5),
# limits=c(11, 13)) +
# scale_y_continuous(name = "Density") +
# theme(plot.title = element_text(size = 14, face = "bold"),
# text = element_text(size = 12)) +
# guides(fill=guide_legend(title=NULL))+
# scale_fill_discrete(labels=c("Control - Public School", "Treatment - Private School"))8.2.2 Efecto tratamiento de matching
Now that we have a matched dataset we can calculate the effect of treatment ( = being in a private school). As before, we can do a simple t-test across the two groups and compare the math skill scores between students in a private school and students in a public school.
# t.test(m1data$math_skills[m1data$school==1], m1data$math_skills[m1data$school==0], paired = T)he t-test shows that there is a statistically significant difference across the two groups (p-value < 0.001). The difference is equal to -0.34, suggesting that students attending a private school report, on average, a math skills score that is 4.33 points lower that students in public schools.
Obtenemos este mismo resultado al realizar la regresión con la dummy de tratamiento:
As before, we can look at results in a graphical form below. The treatment density curve (blue) tend more towards lower scores, marked on the left side of the graph, compared to the control (pink) curve.
# ggplot(m1data, aes(x = m1data$math_skills, fill = m1data$school2)) +
# geom_density(position="identity", alpha=0.6) +
# scale_fill_brewer(palette="Dark2") +
# ggtitle("Density plot of math skill score, by group") +
# scale_x_continuous(name = "Math skill score",
# breaks = seq(0, 100, 20),
# limits=c(0, 100)) +
# scale_y_continuous(name = "Density") +
# theme(plot.title = element_text(size = 14, face = "bold"),
# text = element_text(size = 12)) +
# guides(fill=guide_legend(title=NULL))+
# scale_fill_discrete(labels=c("Control - Public School", "Treatment - Private School"))5. Regresión Discontinua
La última técnica que usaremos el la regresión discontinua. Vamos a hacer un ejemplo muy corto, pero necesitaremos más datos.
Pensemos que para asistir al preuniversitario, los estudiantes tenían que dar una prueba diagnóstico y se invitaron a participar del preuniversitario que los estudiantes solo que sacaron hasta 500 puntos.
IV. Presentación de los resultados
El ultimo elemento que vamos a preparar son un conjunto de tablas que puedan dar cuenta de nuestros resultados.
Para esto usaremos dos paquetes: stargazer & modelsummary